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Definition 
Scattering length 
Atomic mass 
Chemical potential 
Average number of condensate atoms 
Effective two-body interaction strength 
Two-body spin-flip collisional loss rate 
Cylindrical coordinate system 
Trap radial and axial harmonic oscilla- 
tor angular frequencies 



Introduction 



The recent production of trapped atomic vapors of 
87 Rb [1], 7 Li [2], and 23 Na [3] at the phase space densi- 



ties necessary to generate Bose-Einstein condensation 
(BEC) in the ground state of magnetic traps has gener- 
ated significant interest in this area of physics. Previ- 
ously, BEC has only been encountered in the He II phase 
of liquid 4 He, and there has also been accumulating 
evidence for BEC of an exciton gas in cuprous oxide [4] . 
However, it now seems that the fraction of atoms that 
participate in BEC in He II is of the order of 10 %, 
because of the presence of strong atomic interactions 
[5, 6]. In the dilute alkali gas BECs, such interactions 
are much weaker, and experimental techniques have 
been used to produce samples in which the condensate 
fraction is very close to unity. This makes the alkali 
systems most attractive for the study of the essential 
phenomena of BEC, especially since the condensate 
density and temperature, and perhaps even the intrinsic 
atomic interactions [7], are subject to external experi- 
mental control. 

However, the alkali condensates are always produced 
in particle traps, and this endows them with an 
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essentially inhomogeneous character quite distinct from 
the homogensous nature of 4 He II. Detailed understand- 
ing of the properties of these systems requires [8-16] 
that this inhomogeneity be taken into account. The pur- 
pose of this paper is to present simple analytical meth- 
ods which have been found to give results that are in 
good qualitative agreement with those of extensive nu- 
merical calculations based on mean-field theory, i.e., the 
basis-set methods of Refs. [12, 16]. These methods are 
thus expected to be useful to experimentalists who wish 
to explore alternative trap designs, or to modellers who 
wish to analyze large-scale numerical results. 

2. Overview of Dilute Gas Bose-Einstein 
Condensates 

For 87 Rb and 23 Na the scattering lengths, a , are known 
to be positive [3,17,18]. In standard mean- field theory, 
this leads to a stable condensate wave function which is 
larger than the noninteracting ground state of the trap, 
and which can be modeled well by a mean field equation 
[8-16]. For 7 Li, a is known to be negative [2], which 
produces a net effective attractive interaction between 
the particles. A homogeneous condensate with a < is 
predicted to be unstable [19]. However, a spatially con- 
fined condensate has been shown to be metastable [9]; 
and, when the average number of atoms in the conden- 
sate, N , is sufficiently small, the condensate lifetime 
can be quite long. The potential energy of the confined 
system acts as a barrier to the collapse of the condensate 
[9, 11]. Other discussions of the negative scattering 
length case have addressed the general issues of en- 
ergetic stability and the possibility of a transition to a 
denser phase [20-22]. Theoretical predictions of the 
sizes and lifetimes of Bose condensed gases can now be 
compared directly with experiments, as will be dis- 
cussed below. The effect of vortex formation [11-14] 
may also be examined in the present framework. 

In this paper, I examine the general cases of the two 
approximate methods which have been most used to 
study Bose condensed systems, i.e., the Thomas-Fermi 
model [8, 11], and the variational method [11]. While 
the ground state properties of such Bose condensed 
gases have been studied in some detail [8-16], most of 
the methods used have been computationally intensive. 
I show here that approximate methods give good quanti- 
tative results for the ground state properties of these 
gases, and also for vortex solutions of the equations of 
motion. 

3. Basic Results of Mean-Field Theory 

For the low atomic energies and densities achieved in 
the experiments of Refs. [1-3] the basic structure of 



a BEC should be well approximated by mean field 
theory, i.e., the Ginzburg-Pitaevskii-Gross (GPG) 
energy functional [23], described in Eq. (1). In such 
cases the atom-atom interaction is dominated by the 
effects of s -wave collisions. The interaction potential 
may thus be encapsulated in the form V(r, r') = U S 
(r-r'), where U = Air/P'alM. Solutions of the resulting 
many-particle mean-field Schrodinger equation provide 
information about the size, density, chemical potential, 
and lifetime of the condensate as a function of N - The 
evolution of a condensate, after release from the weak 
trap of Ref. [1] for time of flight studies, has been 
described using the time-dependent solutions of the 
Gross-Pitaevskii (GP) equation in Ref. [10]. Here I con- 
sider the condensates formed in the tight trap [1], which 
should be adequately represented by a time-independent 
formalism. The Gross-Pitaevskii-Gross (GPG) energy 
functional [23] for this system may be written 



■/ 



e[i/r] = &r\p*{r) 



■S V2 + ^ 



(r) + -N U \ 



^(r)| 2 hKr), (1) 



where \jj{r) is the common wavefunction for each atom 
in the condensate, and s[ip] is the energy per particle 
in the system. Minimizing s[ip] with respect to varia- 
tions in the wavefunction ip(r), while enforcing the 
normalization condition, Jdr\ip(r)\ 2 = 1, gives the time- 
independent GPG equation [23]. 
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where V tmp (r) is usually an anisotropic harmonic oscil- 
lator potential [1, 2]. Baym and Pethick [11] have iden- 
tified the dimensionless parameter for a cylindrically 
symmetric harmonic system as 



C= (87raN a tmp )' 



(3) 



where a trap = (McOp/A)" 2 , and V tmp (r) = M(co 2 p 2 + 
co 2 z 2 )/2. I use Eq. (1) to model a condensate using a 
variational trial wave function, and in the limit of suffi- 
ciently large £ we will also examine the Thomas-Fermi 
model, using Eq. (2), where the kinetic energy of the 
condensate atoms is neglected. 

An important property, of both current and future 
experimental interest, is the lifetime of the condensate 
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with respect to population loss from the trapped state. 
Such losses derive from spin-flip collisions that force 
the spin-flipped atoms out of the trap. These collisions 
are considered here to be extrinsic to the kinetics of 
BEC formation, and are treated by a simple rate- 
equation approach. The two-body population loss rate 
for a condensate is given by Eq. (20) in Ref. [24]: 



o) = 2a\dr 



R w (Ni 



n 2 (r) + 2n (r)n'(r) + n' 2 (r) 



(4) 



where n (r) and n'ir) are the respective condensate and 
noncondensate number densities, and a is the two-body 
loss rate coefficient associated with spin-flip collisions. 
Since, by the use of forced evaporative cooling, the 
noncondensate atoms may be "cut" away [1], forcing 
n'{r) ~ 0, Eq. [4] reduces to 



R (2 \N ) 



a\Arnl(r) = aNl I dr li//(r)l 4 . (5) 



Finding a solution of the GP equation, or minimum of 
the GPG energy functional, together with the loss-rate 
expression, enables one to estimate the condensate life- 
time in units of a [16, 25]. 

While the ground-state properties of BECs are of 
primary interest, the possibility of vortex formation has 
also been considered [11-14]. A steady-state vortex 
wave function may be written as [26] 



frequency is given by [26] 



h 



IttMR 
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where R is the radius of the system, and £ the radius of 
the vortex core. The critical frequency may be calcu- 
lated directly from the equation AE - L ■ to = 0, where 
AE is the energy gap between ground and vortex states, 
L is the vortex angular momentum, and to is the rota- 
tional angular velocity of the system [26] . In the case of 
a cylindrically symmetric harmonic system, L = N(/n/i , 
and AE = N (s m - s ), where e m is the energy per particle 
in the mth vortex state, and it is a function of N . This 
gives a critical frequency of 



So 



2TTin/i 
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Solving the GPG equation, and finding the energy per 
particle of the condensate, one may calculate the critical 
frequency for formation of such vortices in the system. 

4. Simplified Variational Model 

For a noninteracting cylindrically symmetric har- 
monic system, the eigenfunctions of the Schrodinger 
equation (Eq. (2) with U = 0) are 



<l>n^m(r) = 



e-^ 2l2 H nz (a z z) 



<//(/•) =p»e i 



> '(<•) 



(6) 



where p(r) is the density of the system. Here, S(r), the 
velocity potential, changes by 2 mn along any closed 
loop around the vortex, with m being the vortex winding 
number. For the cylindrically symmetric harmonic traps 
of Refs. [1,2], a wavefunction of the form in Eq. (6), 
with the z-axis as the vortex line, is related to states of 
definite z -component of the total angular momentum. 
The corresponding eigenfunction may be written as 



W') =fn^ z m(r)e 



im<f> 



(7) 



where {n p ,n z ,m} are the radial, axial, and angular 
momentum (z -component) quantum numbers. Equation 
(7) is equivalent to a winding number m vortex state. 
However, for an interacting system, the function f Vzm (r) 
will also be a function of N - As pointed out in Refs. 
[1 1-14], a rotating system should cause the formation of 
a vortex line or ring with critical rotation frequency, j/ crit . 
In a homogeneous rotating system, the critical 
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where a pz = (MQ p J/i) il2 , H n (x) is the n-th order 
Hermite polynomial, and L„ (m) (x) is the m-th associated 
Laguerre polynomial of order n [27]. In the noninteract- 
ing case (i.e., U = 0) Q p = co p and Q z = &> z , however, 
using Eqs. (10) as variational trial wavefunctions for the 
interacting system, with i7 p and O z treated as arbitrary 
variational parameters, we find the GPG energy func- 
tional, 
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where 
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is a dimensionless constant depending on the state being 
considered, normalized so /3ooo = 1, and 
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is important for analysis of vortex states. While Eq. (11) 
is simple to analyze numerically, as I show later, it may 
also be used to give approximate solutions in the limit of 
£ » 1 . While it is easy to examine the perturbative limit, 
I £1 « 1, it is of little experimental interest, since I £ I > 
1 in all experiments to date [1]. 

I minimize Eq. (11) with respect to the variational 
parameters, fl p and fl z , to find, 



2(2t7) 5/2 («, + l/2)nf(o" p 2 (n 2 z -(oj) 
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^iven {2 P = (o p / A , found from dsldfl p = 0, where one has 
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When N is sufficiently large, fl p « co p and O z « co z , so 
one needs only to examine terms of highest-order in co p 
and co z . This is directly equivalent to neglecting the 
kinetic-energy terms in Eq. (11). This approximation 
then gives 
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These expressions may be substituted back into Eq. (11) 
in order to evaluate the energy per particle of the system, 
and one may then calculate the critical frequency for 
vortex formation. For a > 0, the lowest vortex critical 
frequency occurs for the ground state to {n p = 0, n z = 0, 
m = 1 } transition, giving 



von. = t- n » { " p = a " z= a m = 01 - f -V' 5 r 2 , 
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where p 1/e . is the radius of the 1/e population density of 
the ground-state condensate in the z = plane. Unlike 
the homogeneous case of 4 He II, where p ait ~ ln(R/g)/R 2 
from Eq. (8), the trapped condensate has v a - tt ~ MR 2 in 
the large A'o limit, where R = p 1/e is a characteristic size 
of the condensate in the plane of rotation. 

For a < 0, the kinetic energy is never negligible, and 
becomes more important as the condensate population 
approaches the point of collapse. Since the kinetic en- 
ergy may not be neglected in comparison with the self- 
interaction energy, and the potential energy is of roughly 
the same order as the kinetic, the shrinkage of the con- 
densate spatial distribution is not severe enough to make 
the potential energy much smaller than the kinetic en- 
ergy. Hence, Eq. ( 1 4) must be solved numerically. Fig- 
ure 1 shows plots, for the parameters appropriate to the 
experiment of Ref. [2], of the vortex critical rotation 
frequency calculated from numerical solution of Eq. 
(14) for transitions, m = to m = 1 (long-dashed line), 
m = to m = 2 (medium-dashed line), m = to m = 3 
(short-dashed line). Also shown is the frequency for an 
m = to m = 1 transition (solid line), calculated by the 
basis-set method. As can be seen, the lowest critical 
frequency, and therefore the first transition to occur 
when the system is rotated, should be the m = to 
m = °° transition. The maximum stable population for a 
condensate with negative scattering length is given by 
4 = from Eq. (15). Assuming fl z ~ co z , this would give 
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Fig. 1. A plot of the vortex transition frequency, in a trap of frequen- 
cies p p = 163 Hz and v z = 1 17 Hz. The frequencies, calculated by the 
variational method, for m = to m = 1 (long-dashed line), m = to 
m = 2 (medium-dashed line), and m = to m = 3 (short-dashed line) 
transitions to occur. The frequency for a m = to m = 1 transition 
calculated from the basis-set expansion (solid line). 



Fig. 2. A plot of the maximum vortex population (n p = 0, n z = 0), the 
approximate analytic solution Wi ™ 1 = 1453 (m + l)!m!2 2 ™/(2m)! 
(solid line), numerical solution of the minimized Ginzburg-Pitaevskii- 
Gross energy functional (squares), and results from [F. Dalfovo and 
S. Stringari, Phys. Rev. A 53, 2477 (1996)] (circles). 



Equation (21) gives iVSoo ~ 3000 for parameters appro- 
priate to the experiment of Ref . [2] . However, one finds 
N'Sm ~ 1500 by numerical solution of Eq. (14), close to 
the value found in Refs. [9, 12-14]. Figure 2 shows the 
maximum population of a vortex state in the trap of Ref. 
[2], comparing the results of Eq. (22) (solid line), 
numerical solution of Eq. (14) (squares), and results 
from Ref. [13] (circles). As N approaches the maxi- 
mum stable level, many of the properties exhibit 
extreme behavior, i.e., the two-body loss rate has been 
shown to increase rapidly as the condensate population 
approaches the point of collapse [16]. 

The two-body loss rate is given, subtituting Eq. (10) 
into Eq. (5) by 
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Fig. 3. A plot of the two-body loss rate for a 7 Li condensate in a 
cylindrically symmetric trap, v p = 163 Hz and v z = 1 17 Hz. Basis-set 
(solid line), and variational (broken line) calculations. 



This shows, because ^< 1 2 p ) „ zm (A?o)<^ ( m + i)(A?o), that 
vortex states have longer lifetimes than those states 
without, or with less, angular momentum. Figure 3 
shows a plot of the loss rate from a 7 Li condensate in the 
experimental configuration of Ref. [2], the total loss rate 
[including three-body recombination (which was found 
to be negligible except very close to the point of col- 
lapse)] from Ref. [16] (solid line), and the two-body loss 
rate calculated from Eq. (23) (broken line). While the 
variational approach clearly does not show as rapid a 
growth in the loss rate as the critical population, it does 
demonstrate qualitatively the same behavior. 



5. Thomas-Fermi Model 

For a large condensate, t, » 1 , with a > we find 
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meaning the kinetic energy term in Eq. (2) may be 
neglected [8, 11] in comparison to the condensate self- 
interaction energy. This of course is never true for a 
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condensate with a < 0, since the kinetic energy is the 
stabilizing factor in preventing the collapse of the con- 
densate to a more dense phase [9, 11]. For a > 0, how- 
ever, this condition may be satified for a fairly small 
condensate population [8]. Previously, this method has 
been used to study ground state condensates in the large 
N limit [8, 11], but the possibility of studying vortex 
states has only been considered peripherally [12]. While 
Eq. (24) may still be true in the case of a vortex, we may 
not neglect the effect of the angular momentum on the 
spatial distribution of the condensate, since the angular 
momentum imposes a "centrifugal barrier." This 
reduces Eq. (2) to 
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Fig. 4. The number density of a ~2500 atom Rb condensate along 
the radial coordinate, in the z = plane. Basis-set (solid line), varia- 
tional (dotted line), and Thomas-Fermi (broken line) calculations. 



>p{r) = fnf/(r), 



(25) 



where the kinetic energy not related to the angular 
momentum has been neglected, and the angular depen- 
dence, i.e., the phase factor (which does not effect the 
density distribution), of the wavefunction has been 
ignored. For a cylindrical system one has y 2 = m 2 , and 
for a spherical system one has y 2 = 1(1+1), where / is the 
total angular momentum quantum number. This gives a 
Thomas-Fermi wave function of the form 



Ref. [1], where I have assumed that a = 1.3 X 10~ 15 
cm 3 s" 1 [25], comparing the results of the basis-set 
method (solid line), numerical evaluation of Eq. (23) 
(broken line), and the Thomas-Fermi expression in 
Eq. (27) (dotted line). Equation (5) may be integrated to 
give an expression for the lifetime of the condensate 
because of two-body collisions, 
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where V eff (r) = V tra p(r) + h 2 y 2 l 2Mp 2 is the effective trap 
potential. The relationship between N and /x is found by 
the conventional normalization condition, J"drl*//I 2 = 1. 
Figure 4 shows a plot of the number density of a 2000 
atom 87 Rb condensate for parameters appropriate to 
Ref. [1], the wave function calculated by basis-set 
expansion as in Ref. [12, 16] (solid line), Thomas-Fermi 
(brokenline), and variational (dotted line) methods. In 
the simplest case, a harmonic trap with y = 0, many of 
the properties of the condensate may be calculated 
explicitly [8, 11], i.e., the two-body loss rate is given by 
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Fig. 5. A plot of the two-body loss rate for a cylindrically symmetric 
87 Rb trap, v p = 212 Hz and v z = 75 Hz. Basis-set (solid line), varia- 
tional (broken line), and Thomas-Fermi (dotted line) calculations. 
Here a, the two-body loss rate coefficient, is 1.3 X 1CT 15 cm 3 s _1 . 
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where &> 3 rap = u> x u> y u> z = (o 2 p a> z . Figure 5 shows the two- 
body loss rate from the ground state of the trap of 



where T\n,(N ) is the 1/8 population lifetime of a conden- 
sate intially containing Af atoms. This gives a 1/e popu- 
lation lifetime for a 2000 atom 87 Rb condensate in the 
range 32 s to 280 s, for 1.5 X 10' 16 cm 3 s" 1 < a < 1.3 X 
10~ 15 cm 3 s" 1 [25]. This range is very close to the 35 s 
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to 200 s reported in Ref. [12]. The longer maximum 
lifetime of 280 s compared to 200 s occurs because I 
have neglected three-body recombination, while it was 
considered in Ref. [12]. The condensate has a lifetime, 
calculated by the basis-set method, because of spin-flip 
collisions in the range 39 s to 340 s. As we would expect 
the Thomas-Fermi model gives a shorter lifetime than 
the basis-set method, since the condensate is more 
tightly confined (within a finite region of space), c.f. 
Fig. 4. 

I have also calculated the vortex critical frequencies 
for parameters appropriate to the trap of Ref. [1] for 
comparison to the results obtained earlier for the varia- 
tional model, and those presented in Refs. [12-14]. Fig- 
ure 6 shows plots of frequencies for the m = to m = 1 
transition calculated from the large N variational ex- 
pression of Eq. (20) (dotted line), the numerical solution 
of Eq. (14) (broken line), the basis-set expansion (solid 
line), and the Thomas-Fermi model (dot-dashed line). 
For a 1 0000 atom 87 Rb condensate, in the trap of Ref. 
[1], I have found v clil = 26 Hz (basis-set expansion), 20 
Hz (Thomas-Fermi), 18 Hz (numerical solution of Eqs. 
(14, 15), to minimize the GPG energy functional), and 
14 Hz (Eq. (20), the large N approximation), compared 
to the 26 Hz reported in Ref. [13, 14]. While one would 
not expect the Thomas-Fermi approximation to give 
accurate results for relatively small N , as can be seen 
from Figs. (4 and 5), there is reasonable agreement with 
the basis-set calculations for as few as 2000 atoms. 

6. Conclusions 

I have shown that the variational method, using the 
wave functions of Eq. (10), gives very good agreement 
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with the results found by much more computationally 
intensive methods [8-10,12-14,16]. For a > 0, it is clear 
that the variational and Thomas-Fermi methods may be 
used to examine the full range of condensate popula- 
tions, at least when ground states and vortex states are 
being considered, to provide good estimates of the con- 
densate properties, i.e., lower bounds on lifetimes and 
peak densities. The Thomas-Fermi approximation also 
has the advantage of being analytically solvable for the 
most common cases of experimental interest, i.e., har- 
monic traps as in Refs. [1, 2]. It may also be solved 
numerically with little difficultly for traps of more arbi- 
trary geometries [3]. For a < 0, the variational methods 
do not give as good agreement to the numerical calcula- 
tions [13,14,16], especially as the critical population 
levels are approached. While qualitative behavior is re- 
tained, i.e., maximum population levels are in good 
agreement, other properties of the condensate, e.g., two- 
body collisional loss rates, do not show the same 
extreme behavior near the critical population levels. 

The critical frequency for vortex formation in a large 
condensate of repulsive atoms was found to scale as 
Vcrit ~ MR 2 in the limit of large N , where R is the 
characteristic radius of the condensate in the plane of 
rotation, compared to a homogeneous system in which 
I'crit ~ \n(R/g)/R 2 . The formation of vortices, and scal- 
ing of the vortex critical frequency, should be examined 
experimentally in the near future. While the approxi- 
mate methods used here to calculate the vortex critical 
frequencies give reasonable results in comparison to the 
methods of Refs. [12-14], the scaling behaviour at large 
N appears to be very similar for all the methods. In the 
case of negative scattering length, the formation of vor- 
tices has been considered as a mechanism for forming 
condensates containing greater numbers of atoms 
[13, 14]: an m = 1 vortex may contain —4000 atoms, 
compared to ~ 1300 atoms for the ground state in the 
trap of Ref. [2]. I have shown, however, that the forma- 
tion frequency for such vortices is lowest for the m = °° 
state, so the actual formation of such structures in exper- 
iments must be studied more fully. Indeed a time-depen- 
dent calculation may be required to examine the conden- 
sate response to rotation near the critical vortex 
formation frequency. It is possible that vortices may be 
formed during the evaporative cooling cycle, as high 
energy atoms are "cut" away, and low energy atoms are 
ejected because of the rapid increase in loss processes 
[16,20-22] for lower energy states. 



Fig. 6. A plot of the vortex transition frequency, for m = to m = 1 
with a positive scattering length: the variational value calculated 
numerically (broken line), the variational expression for large N 
(dotted line), the Thomas-Fermi approximation (dot-dashed line), and 
the basis-set expansion method (solid line). 
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